Gapped phase in AA stacked bilayer graphene. 
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AA-stacked bilayer graphene supports Fermi circles in its bonding and antibonding bands which 
coincide exactly, leading to symmetry-breaking in the presence of electron-electron interactions. 
We analyze a continuum model of this system in the Hartree-Fock approximation, using a self- 
consistently screened interaction that accounts for the gap in the spectrum in the broken symmetry 
state. The order parameter in the groundstate is shown to be of the Ising type, involving transfer 
of charge between the layers in opposite directions for different sublattices. We analyze the Ising 
phase transition for the system, and argue that it continuously evolves into a Kosterlitz-Thouless 
transition in the limit of vanishing interlayer separation d. The transition temperature is shown to 
depend only on the effective spin stiffness of the system even for d > 0, and an estimate its value 
suggests the transition temperature is of order a few degrees Kelvin. 

PACS numbers: 61.46.-w, 73.22.-f, 73.63.-b 



INTRODUCTION 



Graphene is a two-dimensional triangular lattice of car- 
bon atoms with two sublattice atoms, A and B, in each 
unit cell, forming a honeycomb structure. Its electronic 
low energy properties are governed by a massless Dirac 
Hamiltonian, and near the neutrality point the quasi- 
particle energies disperse linearly with a speed vf OP- 
Coulomb interactions in this system can be character- 
ized by an effective fine-structure constant j3 — e 2 /ehvF, 
where e is a dielectric constant due to the substrate upon 
which graphene is deposited. For moderate values of 
P, electron-electron interactions renormalize the speed 
of the carriers near the Dirac points, but do not alter the 
semi- metallic character of the spectrum [21 13] . For larger 
values of (3 (/? ~ 1), theoretical studies indicate that the 
Coulomb interaction becomes important [3- 6 J , and there 
is the possibility of a gap opening in the spectrum due 
to dynamical symmetry breaking. At present, there is no 
evidence for such a gap under experimentally realizable 
circumstances. 

The situation becomes dramatically different when one 
considers two layer graphene systems [7j . Under most cir- 
cumstances, graphene bilayers are found to have Bernal 
stacking (as in their parent material graphite). Labeling 
the layers as left (L) and right (R), the Bernal structure 
has atoms of one sublattice in the left layer (AL) adjacent 
to atoms of the other sublattice in the right layer (BR) so 
that these are tunnel- coupled, while the other two atoms 
in the unit cell [left layer, B sublattice (BL) and right 
layer, A sublattice (AR)] have no interlayer coupling. In 
this configuration the energy varies quadratically about 
the Dirac points. The large density of states associated 
with this leads to broken-symmetry groundstates even 
for an infinitesimal electron-electron interaction. Many 
such states have been proposed for the system [51420) . 
and some transport experiments suggest a gap opening 
in the spectrum [2D - I24] as is expected in many of these 
scenarios. 
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FIG. 1. (Color online) Energy spectrum for AA stacked bi- 
layer graphene. Solid lines correspond to the Hartree-Fock 
energy bands for /3=1.2 near k c . Dashed lines are single elec- 
tron energies in the absence of interactions. Inset: Nonin- 
teracting single electrons energies over a broader momentum 
range. Solid (dasehd) lines represent bonding (antibonding) 
states. 



Other stackings for bilayer graphene are possible. Re- 
cently, AA stacked bilayers have been identified experi- 
mentally [25] , in which all atoms of one layer are directly 
adjacent, and tunnel coupled, to the equivalent atoms in 
the other layer. This type of stacking also occurs locally 
in large regions of small angle twisted bilayers [5S[ [57] • 
The focus of our study will be on the electronic states of 
such AA stacked graphene bilayers. 

When interactions are ignored, the system may be con- 
veniently represented in terms of bonding and antibond- 
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ing states, each separately supporting its own Dirac spec- 
trum, with the Dirac points separated in energy by 2t\, 
with t\ the interlayer hopping parameter. At half filling 
the Fermi surface coincides with the circle in momen- 
tum space where these spectra cross. When interactions 
are introduced, this spectrum becomes unstable to the 
opening of a gap due to the perfect nesting of the Fermi 
surfaces in the two bands. This situation is analogous 
to what occurs in double layer graphene separated by a 
thin insulating barrier, in which tunneling between lay- 
ers is suppressed but interlayer interactions are not. It 
is widely believed that a gap opens in the spectrum in 
this system, with the groundstatc forming an interlayer 
exciton condensate [28H3T] . In the limit of vanishing in- 
terlayer separation d, the two systems indeed are isomor- 
phic. 

When d is non- vanishing, the double layer graphene 
system and the AA stacked bilayer differ in an important, 
qualitative way. The double layer Hamiltonian has U(l) 
symmetry, which is broken in the groundstate, whereas, 
as explained below, the AA bilayer system has only a Z2 
(Ising) symmetry which is also broken in the groundstate. 
In what follows we compute the bandstructure of this sys- 
tem in the Hartree-Fock approximation, demonstrating a 
gap opening associated with the broken symmetry. This 
gap serves as an estimate of the mean-field transition 
temperature. To take into account the subtle effects of 
screening, we include a static dielectric constant in the 
effective electron-electron interaction, which is computed 
self-consistently in the presence of the gap in the energy 
spectrum [31]. For typical values of the system param- 
eters, the resulting gap is of order 5meV, suggesting a 
mean-field ordering temperature of 50K. In the broken 
symmetry state, we find a staggered charge density dis- 
tribution such that charge is transferred between layers, 
but in opposite directions for each of the sublattices. Al- 
though our calculations do not include spin explicitly, it 
is likely that this charge density wave ordering would be 
compensated by the two spin species, yielding antiferro- 
magnetic spin ordering in the system [32] . 

In two dimensions the mean-field transition temper- 
ature greatly overestimates the true disordering T c by 
missing the low energy collective excitations, as well as 
crucial topological excitations. In the limit d 0, where 
the system has U(l) symmetry, the latter are vortices, 
and the disordering transition should be in the Kosterlitz- 
Thoulcss (KT) universality class [33]. In this case the 
transition temperature is controlled by the phase stiffness 
of the U(l) degree of freedom. This stiffness is challeng- 
ing to compute accurately because it depends exponen- 
tially on the effective interaction strength, so that differ- 
ent assumptions about screening lead to vastly different 
estimates of the critical temperature [55] [30] [3TJ [3H [35] . 
For d > 0, the disordering transition falls into the Ising 
universality class, and may be understood as being driven 
by a proliferation of domain walls [33] - Despite the dif- 
ference from the KT transition, we argue below that the 
transition is still controlled by the effective phase stiff- 



ness, and that in principle the Ising transition would con- 
tinuously go into the KT transition if the layer spacing 
d could be continuously tuned to zero. Combining this 
analysis with a calculation of the phase stiffness using 
the self-consistently screened interaction, we arrive at an 
estimate of T c of order 3K for the AA stacked graphene 
bilayer. 

In principle, the Ising transition that we analyze in 
this work should present itself as singularities in ther- 
modynamic quantities of the system. These would be 
considerably less subtle than what is expected for the 
KT transition associated with the related double layer 
system. More directly, the symmetry-breaking in the 
groundstate should be visible at low temperature in STM 
experiments with sufficient resolution to distinguish the 
two sublattices. 

This article is organized as follows. In Section II, we 
describe the model used to describe the system, a sin- 
gle Dirac cone treated in a continuum approximation. 
In Section III we describe the broken symmetry state 
within the Hartree-Fock approximation, including the 
self-consistent approach to screening of the interactions. 
Section IV describes a simple continuum model of the 
energetics of the order parameter, how this leads to an 
Ising transition controlled by the phase stiffness in the 
energy functional, and the resulting estimate of the dis- 
ordering temperature. We conclude in Section V with a 
summary of our results. An Appendix provides details 
of how screening is incorporated into our calculation. 

II. MODEL 

We consider a single valley model for electrons in 
graphene in the continuum limit. Near the Dirac point 
the non-interacting Hamiltonian takes the form [36] 

H = hv F ^ (k x - iky )c\ ik c B ik - h cl Lk c am + h.c., 

V,i=L.R k,o=A,B 

(1) 

where the operator c aik creates an electron on sublattice 
a in layer i with momentum k. The interlayer hopping 
matrix element is t\, which is estimated to have values 
of order ~ 0.12t [7], where t is the nearest neighbor in- 
tralayer hopping parameter. In the continuum limit, t en- 
ters only through the electron speed HvF = y ^ta, where a 
is the monolayer graphene lattice parameter. Valley and 
spin degrees of freedom are not explicitly treated. The 
eigenvalues of Hq are e(fc)±ii, and represent bonding (-) 
and antibonding (+) linear combinations of left and right 
layer wave functions. The monolayer graphene dispersion 
is e(k) — shvpk, where s = —1 indicates an electron-like 
band and and s = +1 a hole-like band. 

The inset of Fig{T] illustrates the one-electron band 
structure as a function of the momentum near a Dirac 
point. The result is highly analogous to what one finds 
in the biased double-layer case 2* 301. As in that case, 
when the system is neutral the Fermi surface corresponds 
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with the circle of radius k c — ti/hvp at which the hole- 
like antibonding band spectrum crosses the electron-like 
bonding band spectrum. Because of this perfect nesting, 
the system is unstable to the formation of a particle-hole 
pair condensate 37, 38] when interactions are included. 

We model the interaction part of the Hamiltonian as 



< q E/ E/ E/ ' C o='' 



25 



k'+q 



a,/3 i,j q,k,k' 



(2) 

where 5* is the sample area, V LL (q) = V RR (q) and 
V LR (q) = V RL (q) are Fourier transforms of a screened 
intra- and inter-layer electron-electron interaction, re- 
spectively, and : ... : indicates normal ordering. Eq[2] 
assumes that the interactions preserve the sublattice and 
valley indices [3J. 



100 



80 



§ 60 

"-I — i 

o 

B 40 



20 



\ P=1 .2 
\ Gap=0.005t 
\ 

\ 

\ 

s 



V LL (q) 



V LR (q) 

Bare Coulomb 



0.0 0.2 0.4 0.6 0.8 1.0 



Because of the perfect nesting of the Fermi surfaces in 
the bonding/antibonding bands, we expect that electron- 
electron interactions drive an instability involving the 
mixing of these bands near the Fermi circle at k c . Within 
the Hartree-Fock approximation this is expressed as a 
condensation of electron-hole pairs. Any estimate of the 
electron-hole pairing amplitude depends sensitively on 
the level of screening included in the model. Because of 
the long range character of the Coulomb interaction, an 
unscreened electron-electron interaction results in very 
large gaps in the spectrum, and in a high critical tem- 
perature |28l I29] , To go beyond this, one may incorpo- 
rate screening via a dielectric function appropriate for 
metallic graphene, completely eliminating the effective 
long-range nature of the interaction. This results in both 
a very small gap and a very low critical temperature 
[501 SO]. In this work we adopt a model of static 
screening [UJ which allows for its suppression at large dis- 
tances, due to the appearance of the gap in the spectrum 
due to electron-hole condensation. The gap and dielec- 
tric screening functions are computed in a self-consistent 
way [3TJ[35 ( 35 • hi the Appendix we provide some details 
of the polarization function matrix and how this leads to 
our model of the screened potential. 

In Figj2] we plot a typical result for the inter- and 
intra-layer interactions for a condensate with assumed 
gap Eg = 0.005i. At large q the interactions are highly 
reduced from the bare Coulomb interaction. However, 
for values of q below ~ E g /hvp, V (q) and V LR (q) are 
not fully screened, and one finds a Coulomb- like 1/q di- 
vergence. This model thus produces an intermediate be- 
havior between approaches using unscreened and metal- 
lically screened interactions. A complication of this ap- 
proach is that it requires knowledge of the single-particle 
band structure, so that the screened interactions must be 
computed self-consistently with states and energies found 
in the Hartree-Fock approximation. We next turn to a 
discussion of this calculation. 



FIG. 2. (Color online)Inter- and intra-layer interaction as 
a function of q, for a condensate with /3 = 1.2 and gap 
E g =0.005t. For comparison the bare Coulomb interaction 
is also shown. 



III. 



HARTREE-FOCK APPROXIMATION 



In the Hartree-Fock (HF) approximation the electron- 
electron interaction is replaced in the Hamiltonian by 



ttHF _ 

ti„„ — — 



-5V 



k c,9 jk F lJ (|k-k'|) < c f k ,c mk ' > 



(3) 

The energy of the groundstate can be lowered by opening 
a gap at the crossing between the electron-like band of 
the bonding energy states and the hole-like band of the 
antibonding spectrum. To do this we consider broken 
symmetry states characterized by the circularly symmet- 
ric order parameters 



A(k) 



iB(k) = < c^^CALk >=< c ] BRk c B Lk >* and 
Q(k) = < c ARV .c A Rk >=< c BLk c B LU > 



1- 



J 



~< C BBk C BRk >= !-< c ALk°ALk > (4) 

These forms are valid provided one assumes the electron- 
like antibonding band is completely empty, and the hole- 
like bonding band is completely full. These assump- 
tions are sensible because neither of these bands approach 
the Fermi energy in the non-interacting state. The self- 
energies associated with these parameters are 

1 



A L R R (k) +l Af R (k}= 
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(Ik- 
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Rewriting Eq|3] in terms of these quantities, the re- 
sulting Hartree-Fock Hamiltonian may be diagonalized, 
yielding eigenenergies ±e±(fc) of the form 

e± (k)=^/(hv F k t {h+At R (kW + (Af R {k) + A LL (k)) 2 . 

(6) 



In terms of the self-energies and eigenenergies, at finite 
temperature the order parameters take the forms 



A(k)+iB(k) = - 



-hv F k+ti- Ak R (k)-HA^ R (k) , (e(k) 

■ tanh 1 



\2k B T 



r 



hvph+t-y- A^ R {k)+iAf R (k) 



tanh 



2k B T / 
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Q(A)= 
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FIG. 3. (Color online) Order parameters as a function of k for 
/3=1.2. Inset: Schematic plot of the direction and magnitude 
of the isospin on atoms A and B as a function of momentum. 

In Fig[T]and Fig(3]we plot the band structure and the 
order parameters obtained from self-consistently solving 
Eqs. [4][8]for /3=1.2. As expected the coherence between 
the bonding and antibonding bands opens a gap near k c . 
Moreover, the minimum gap occurs at larger momentum 
than k c because the self-energy A^ fl acts as an extra 
contribution to the hopping parameter connecting the 
layers. 

The inset of Fig[2] illustrates the value of the gap as 
a function of the strength of the Coulomb interaction j3. 
In this example we adopted a hopping matrix element 
t ~ 2.8eV . A moderate value of j3 ~ 0.6 corresponds to a 
background dielectric constant of ~ 3.7, which could for 
example be provided by a Si02 substrate. At this value 
of P one hnds E g ~ b.QmeV, suggesting a mean-field 



transition temperature of several tens of degrees Kelvin. 
However, as we discuss below, fluctuations will consider- 
ably reduce the critical temperature. 

The ordering in this system arises from spontaneous 
coherence between independent bands. For momenta 
k -C k c , the two occupied bands have a bonding char- 
acter and A(k) rj 1/2. For values of k ^> k c , the hole- 
like antibonding band and the electron-like bonding band 
are occupied, and A(k) m 0. In these two limits Q(k) is 
close to its non-interacting value, Q(k) ~ 1/2 because 
the large energy differences between states in different 
bands allows for little admixture between them in the 
HF state. For values of k ~ k c , the similarity in en- 
ergy of bonding and antibonding states allows them to 
admix and thereby lower the energy of the system. It is 
in this range that Q(k) significantly differs from its non- 
interacting value, so that we can use this quantity as a 
measure of the anomalous order in the HF state. 

A physical picture of the ordering can be found in 
the language of quantum ferromagnetism. We define an 
isospin degree of freedom for each of the sublattices, in 
which spin "up" corresponds to an electron in the left 
layer, and spin "down" indicates an electron residing in 
the right layer. Assuming we consider only states in 
which the hole-like bonding band is completely full, and 
the electron-like anti-bonding band is completely empty, 
one may easily confirm that the isospin expectation val- 
ues in the HF state can be written in terms of spin vectors 
for each sublattice, satisfying 



S x ,A(k) 

SyAk) = 

S z ,A(k) = 



S x , B {k) 
-S y , B {k) 



2A(k) 
2B(k) 
2Q(k) 



1. 



(9) 



Note that, because of the background filled hole-like 
bonding band, these vectors are not constant in mag- 
nitude. However they can be written in terms of a vector 
field of fixed magnitude c?(fc), with constant vectors in the 
x direction subtracted (added) for the A (B) sublattice. 

In the non-interacting state, one finds S x for both sub- 
lattices takes on its maximal value for k < k c , and Sf ytZ ) 
both vanish. As k increases through k c , there is a sudden 
change in which S collapses to zero. The introduction of 
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nature of the Coulomb interaction is included. 
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FIG. 4. (Color online) Difference in energy per unit area, as 
function of the strength of the electron-electron interaction, 
between the Ising and the U(l) phases. See text. Inset: Gap 
as a function of /3. 



the non- vanishing Q order parameter allows S to evolve 
continuously from its small k value to zero at large k, 
lowering the isospin exchange energy and thus the total 
energy of the state. This behavior is illustrated in Fig. [3] 
In this evolution, a rotates through the y direction with 
no component in the z direction, because tilting in the 
latter direction amounts to a charge imbalance between 
the layers, adding a positive charging (Hartree) energy 
to the total groundstate energy. Note that in the limit 
d — >• 0, this latter energy cost vanishes, so that the rota- 
tion can be through any direction in the y — z plane with 
no additional energy cost, and the groundstate has a bro- 
ken C/(l) symmetry. For non- vanishing d 1 the only choice 
that does not affect the energy of the state is whether the 
rotation occurs through the positive or negative y direc- 
tion. One thus obtains a broken Z 2 symmetry in the 
groundstate. 

Physically, this broken symmetry manifests itself as a 
transfer of charge between the layers, in opposite direc- 
tions for the two sublattices. A related ordering for this 
system has recently been suggested in a tight-binding 
Hubbard model [35] • Were spin to be included in our 
model, we expect that the spontaneous orderings for op- 
posing spin directions would be staggered, yielding the 
essentially the same antiferromagnetic ordering in our 
system that is found in the tight-binding Hubbard model. 
We note that the very short-range interaction employed 
in a Hubbard model implies strong screening, which is 
not consistent with the gap opening in the spectrum. 
Our calculation strongly suggests that the antiferromag- 
netic ordering will still be present when the long-range 



IV. FINITE TEMPERATURE ISING 
TRANSITION 

In AA stacked bilayer graphene, the distance between 
layers is relatively small (d=\.2a, with a the graphene 
lattice constant.) Moreover, screening is very efficient 
in this system at large q. In the range of wavevectors 
relevant for the formation of the condensate (fc ~ k c ), 
these two effects lead to intra- and inter-layer electron- 
electron interactions which are very similar, as illustrated 
in Figj2j Because of this, the difference in energy be- 
tween a phase with U(l) symmetry obtained by tak- 
ing V LL (q)=V LR (q), and one with broken Ising symme- 
try, from a Hamiltonian in which the difference between 
the intra- and inter-layer interactions is retained, is very 
small. This is illustrated in Fig(4j It is then interesting 
to consider how the KT transition expected in the d — > 
limit is related to the Ising transition expected for the 
d > case. In this section we argue that the transitions 
are continuously connected, based on a scaling argument, 
and that both transition temperatures are set by the ef- 
fective phase stiffness of the system. 

To address this issue, we consider an effective low en- 
ergy functional of the form 



E 



dr 



-|V(9(r))| 2 - ftcos(20(r)) 



(10) 



where 9 is a phase angle for the electron-hole condensate, 
which can be interpreted as the angle in the y — z plane 
made the vectors S a described in the last section, h is 
the charging energy cost for d > 0, and p s is the phase 
stiffness of the U(l) degree of freedom. For h = the 
quantity p s may be estimated by |42j 



p. CO 



ii 



-k B T 



E 



e oo [(7rfcBT(2n + l))2+ J B g)0 2 (T)]3/2 

where E g0 is the energy gap when d = 0, which is some- 
what smaller than E g . 



Eq. 10 may be placed on a lattice to allow for the pos- 
sibility of vortex excitations. For h = the Hamiltonian 
then has the form Hxy = ~JJ2<ij> cos [#i ~ where 
0i is a phase angle at site i, which is assumed to be on 
a square lattice, the sum is over nearest neighbors sites, 
and J ~ p s . This is just the XY model, and it supports 
a KT transition at T KT = 0.89 J Eq. [To]for h = 
and Hxy share the same U(l) symmetry, and provided 
one allows for vortex configurations in the former, both 
models support a KT transition. 

This observation is dramatically changed when h > 0, 
which introduces a perturbation in the Z7" (1) theory that 
is relevant for any temperature below Tkt [44] - The 
usual interpretation of this behavior is that the KT tran- 
sition is eliminated 45, 46J, and the only thermodynamic 
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phase transition remaining is of the Ising type. In apply- 
ing Eq. 10 to the AA bilayer graphene problem, the 
smallness of the effective h, which vanishes in the d — » 
limit, naively suggests a very small T c associated with 
this Ising transition. With Tkt ~ Ps> if this is so then it 
is unclear what becomes of the KT transition as soon as 
h is raised from zero. 

To address these questions we first consider how to 
construct an Ising model that appropriately describes 
Eq. 10 as the Ising transition is approached. This 



the well-known form 



can be accomplished by noting that the Ising transition 
may be interpreted as a proliferation of domain walls 
[33], so that one should match the energy cost of over- 
turned neighboring Ising spins er z in the Ising Hamilto- 
nian H Ising = -K J2<ij> <^z{i)^z(j) with the energy of 
a domain wall (DW) in Eq. [ToJ Moreover, the nearest 
neighbor distance in Hj S i ng should be taken as the DW 
width of Eq. [10] The DW configuration of Eq. [I0]may be 
explicitly computed with standard methods [33] , yielding 
the form 



9(x) — ±2 arctan{exp[x/£]}, 

with an energy per unit length Sdw = 4^/p s /i and DW 
width £ = y p s /h/2. Identifying this last quantity with 
the lattice constant of Hj sing , we arrive at the esti- 
mate 2K = £dw£, — 2p s . Using the two-dimensional 
DW proliferation transition temperature |33] ksTpw — 
2K/ log(2.63), we arrive at an estimate for the Ising tran- 
sition in the bilayer system of fcsTj S j ng = 2p s / log(2.63). 

Remarkably, h does not explicitly enter into this re- 
sult, so that the transition temperature does not vanish 
as h — > 0. Physically, this can be understood by recog- 
nizing that for small h, the unit cell area £ 2 of Hj sing 
encompasses an increasingly large physical area of the 
system as h — > 0, such that the energy to overturn 8 over 
this length scale remains finite even as the microscopic 
energy cost to do so vanishes. It is interesting to note 
that our estimate of ksTising, like ksTxT, depends only 
on p s , and in fact exceeds the standard KT transition 
temperature, which in the simplest estimate |33j is given 
by k B T KT — ~p s . We are then left with the question of 
how the Ising transition vanishes and the KT transition 
emerges as h — > 0. 

The key to answering this question is to note that the 
two transitions involve the proliferation of topologically 
distinct defects, so that relating them requires a model 
that has both of them. A simple way to do this is to 
assume that the perturbation h may be ignored below 
the scale £, so that p s is renormalized by vortices in just 
the way it would in the absence of the U(l) symmetry- 
breaking term; above this scale, we assume the vortices 
are linearly confined into vortex-antivortex pairs [47] , and 
can be ignored at such long length scales [4~5j . 

To carry this out, we use the renormalization group 
(RG) flow equations for the KT transition, which have 



dC 
dy 



d£ 



= [2 - Trp s ]y 



(12) 
(13) 



where y is the vortex fugacity, assumed to be small, 
p s = ps/ksT, and I characterizes the multiplicative 
rescaling factor of the shortest length scale in the effec- 
tive Hamiltonian, a — > ae l . We assume the ultimate scale 
for this microscopic length scale is £, which is connected 
to the phase stiffness by p s = (/i//cbT)£ 2 . 



log(2.63) 




FIG. 5. (Color online) Schematic RG flows of p s as a function 
of scale size f. Flow stops when p B (£) matches h£ 2 /ksT. This 
value sets the exchange energy for an effective Ising model. 
Dotted (blue) line indicates critical value for Ising transition. 
Dashed (red) line indicates critical value for KT transition. 

Fig. [5] illustrates this relation between (the renormal- 
ized) p s and £, along with the RG flows for different initial 
values of p s . For simplicity we ignore the renormaliza- 
tions of y and h; including these would not alter the 
qualitative physics. The initial value of p s (from Eq. 10 ) 



sets the value on the furthest left of the flow, and as can 
be seen in the figure, p s is always renormalized down- 
ward. Following the logic described above, the value of 
K used in the effective Ising model should be taken as 
K = ksTp s (£,), with £ satisfying the self-consistency con- 
dition p s (0 = h£ 2 /k B T. 

The KT RG flows have a well-known property, that 
in the limit £ — > 00, p s jumps from 2/-7T to zero as T 
passes from below to above Tkt- We illustrate this crit- 
ical value as a dashed line in Fig. [5] Assuming the Ising 
transition occurs, as remarked above, at a higher transi- 
tion temperature than the KT transition for h = 0, the 
critical p s at which the Ising transition occurs is below 
the critical value for the KT transition. This is illus- 
trated as dotted line in Fig. [5] One can now see how 
h controls the Ising transition. For relatively large h, 
h£, 2 /ksT will curve sharply upward, and p s will undergo 
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only a small downward renormalization, and we expect 
kBTismg ~ i Qg ^2°63) = 2 - 07 /°s, with Ps close to the initial 
value appearing in Eq. [TUj As h decreases, the downward 
renormalization of p s becomes increasingly pronounced, 
lowering the critical temperature for the Ising transition. 
In the limit h — > 0, p s cither flows to a value above 2/w, 
which is above the dotted line (indicating an ordered 
Ising state), or flows to zero, below the dotted line (in- 
dicating a disordered Ising state.) In this way, the Ising 
transition continuously flows into the KT transition in 
the limit h — >■ 0. 

This indicates that, in an h — T phase diagram, the 
transition will be in the Ising universality class for any 
non- vanishing h, but at the endpoint h = 0, the transi- 
tion will be in the KT universality class. 

In Fig. [6] wc illustrate an estimate of the critical Ising 
temperature for the AA bilayer graphene system when 
vortices are ignored. In this situation the mean-field 
Ps(T) may be computed from Eq. the Ising tran- 
sition will then approximately occur at the temperature 
where this coincides with p s (T) = k B T/2.07. Fig. [6^a) 
illustrates this matching condition, and [G^b) illustrates 
the resulting transition temperature as a function of the 
electron-electron interaction strength /?. For f3 w 0.6 we 
arrive at an estimate of T c « 10~ 4 £//cs ~ 3K. 




0.0 0.2 0.4 0.6 
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V. SUMMARY AND DISCUSSION 

In this work we have considered broken symmetry 
states of AA-stacked bilayer graphene. Due to the inter- 
layer tunneling this system has two Dirac cones symmet- 
rically displaced above and below zero energy, so that at 
zero doping and in the absence of interactions, there is a 
band crossing precisely at the Fermi circle. This perfect 
nesting leads to an instability in the presence of arbi- 
trarily small interaction strength, in which a gap opens 
around the Fermi energy. This gap arises from sponta- 
neous coherence between the bonding and antibonding 
bands, in a manner analogous to what happens in dou- 
ble layer graphene systems with an interlayer bias. This 
spontaneous coherence leads to a transfer of charge from 
one layer to the other, in opposite directions for the two 
sublattices so there is no net charge imbalance between 
layers. 

Unlike the double layer case, the difference between 
inter- and intra-layer interactions in the AA-stacked bi- 
layer leads to an Ising-like order parameter, rather than 
one with an XY character. We argued that the Ising 
transition in the former case and KT transition in the 
latter are continuously connected as the layer spacing 
d vanishes (eliminating the difference between the two 
types of interactions.) Both transitions are controlled by 
the effective (iso-)spin stiffness of the system. A calcu- 
lation of this, including screening which self-consistently 
incorporates the effect of the gap, yields a critical tem- 
perature estimate of a few degrees Kelvin. 

Our analysis neglects the presence of multiple Dirac 
points of the single layer system, which will be present 
due to valley and spin degeneracy. Incorporation of this 
effect is likely to yield multiple copies of the order param- 
eter in the system; the spin degree of freedom in particu- 
lar should allow the charging pattern to be neutralized in 
by having the order parameter for opposite spins pointing 
in opposite directions. The resulting groundstate would 
have antiferromagnetic order [32] . 

In principle, the staggered groundstate ordering could 
be detected by spin-polarized tunneling into one of the 
two sides of the structure. It is also interesting to note 
that thermodynamic measurements could, in principle, 
show signs of the expected Ising transition, for example 
in singular behavior of the heat capacity. Such thermo- 
dynamic singularities are considerably more accessible in 
this system than in the double layer system because dis- 
ordering in the latter occurs through a KT transition, 
which is considerably more subtle that the Ising transi- 
tion. 



FIG. 6. (Color online) (a) Graphical solution of equation 
kaTising « 2.07p s for the critical temperature. Continuous 
line corresponds to p s (T) and the red dashed line is a straight 
line of slope 2.07. (b) Values of the critical temperature as 
function of the strength of the electron electron interaction. 
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VII. APPENDIX: DIELECTRIC CONSTANT AND SCREENED INTERACTION. 



We define the polarizabilitics 



9s g v 
S 



E 

k,s,s 



n s ,k - »s',k' 



i^(k,k') 



(14) 



where k'=k' + q, g s = 2 and g v = 2 are the spin and valley degeneracies, S 1 is the sample area and e s ^ are the 
eigenvalues of the band s with wavevector k. The indices a and (3 are the layer indices. The couplings F have the 
form 



F fl Q ;?(k,k') =< S ,k\P a \s'y >< s',k'\Pp\s,k >, 



where P a is a layer projection operator, 

Pl = 



/l 000 

10 



\ 



and Pr 



/ \ 



10 

\ 1 / 



(15) 



(16) 



These operators are written in the basis LA,LB,RA,RB. Eq. 14 depends on both the eigenvalues and eigenvectors 
of the HF Hamiltonian, which in turn depends on the order parameters characterizing the condensate. By symmetry 
= Hr,r and Xl L ,R = Kr,l- 
Using these polarizabilities we obtain effective intra- and inter-layer interactions in the RPA 



V LL (q) V LR (q) 
V LR (q) V LL {q) 



v(q) 



1 
1 



v(q) 



1 



,-qd 



where d is the distance between layers. From this equation we finally obtain 

I + v(q)Tl LL (q)(l - e~ 2 i d ) 



n L L(<z) Il LR (q) 

iW<z) n LL ( g ) 



-qd 



l 

3 -qd 1 



(17) 



V L L{q) = v{q) 
V L R(q) = v(q) 



(1 + v(q)IL LL (q) + v{q)Tl LR (q)e-<i d y - (v(q)U LR (q) + v(q)e'i d U LL (q)y 

e~o d + v{q)Tl L R{q)(e~ 2 i d ~l) 

(1 + v{q)U LL (q) + v{q)U LR {q)e-idy - (v(q)U LR (q) + v{q) e -i d Il LL {q)) 2 ' 



(18) 
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